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Abstract 

Thermal  diffusion  in  stratified  media  is  important  in  many  applications  such  as  electronic 
packaging,  optical  disk  recording,  and  nondestructive  evaluation.  A straightforward 
algorithm  is  presented  for  the  solution  of  the  one  dimensional  thermal  diffusion  equation  in  a 
stratified  medium  containing  a modulated  heating  source.  In  this  context,  media  with  thermal 
conductivity  (or  diffusivity)  gradients  normal  to  the  specimen  surfaces  can  be  considered 
stratified  media.  Interface  thermal  resistance  is  easily  incorporated  into  the  calculation.  In 
the  case  of  heat  sources  varying  in  the  x-y  plane  parallel  to  the  layers,  the  procedure  can  be 
used  to  find  the  x-y  Fourier  transform  of  the  temperature.  In  the  case  of  a multilayer  stack 
of  rectangular  cross  section,  the  procedure  can  be  used  to  fiind  the  coefficients  of  the  x-y 
Fourier  series  expansion  of  the  temperature.  A sample  calculation  for  a delaminated 
diamond  film  on  a WC  substrate  is  presented. 

Key  words:  diamond  film,  heat  diffusion,  multilayer,  stratified  medium,  thermal  diffusion, 
thermal  waves 


1 


T ■ 


S;4’ 

y mk 


A Hi  iiutTAWOd  wwariwa  oAMiHsmaMT  ^ morTum^liftem* 
<i’.>M'(  io«  viMPfAan  ^A>aN4<?s»<J  A ifrr»  " 

^ " ■'  ^ 'ffV^  - ^V’  i?!' 

ne^feAf  ^■■-  - ' ■ ^ 

iim  to’ 

?maM  .(jTf»te»ittoje  ,_#'*■•%  ' 


f*r 


■M 


'li , u til  yo  9^3 , 

^ >^*^irmami:ifOo  9ffy  |^‘o/  ,0hom 


■5'  ’ ^ 


■« ' 'tfeomtfir  ■ 


,'  1- 


:m 


m. 


, Sj  v,!! 


M' 


1^' 


jj^4a’  'I* 


:M' 


5, 


re 


.L 


ALGORITHM  FOR  SOLUTIONS  OF  THE  THERMAL  DIFFUSION  EQUATION  IN  A 
STRATIFIED  MEDIUM  WITH  A MODULATED  HEATING  SOURCE 


1.  Introduction 

The  solutions  of  the  thermal  diffusion  equation  in  a stratified  (multilayer)  medium  are 
important  for  many  modem  applications  such  as  thermal  diffusivity  measurements  ^ thermal 
wave  nondestmctive  probes^,  heat  conduction  in  optical  disk  recording^,  and  heat  dissipation 
in  electronic  multilayer  stmctures"^.  Several  analyses  have  been  made  of  ac^  or  dc  steady 
state  heat  flow  in  general  multilayer  systems^'^,  however,  they  have  included  either  no  heat 
source^  or  a heat  source  at  an  outer  surface^ 

This  paper  presents  a simple  algorithm  to  solve  the  one  dimensional  heat  diffusion 
equation  for  a general  multilayer  stack  containing  a modulated  (ac)  heating  source  located  at 
an  arbitrary  plane  within  the  stack.  The  algorithm  used  involves  a simple  matrix  formalism 
that  is  analogous  to  that  given  in  Carslaw  and  Yaeger^.  Thermal  resistance  at  each  interface 
is  easily  incorporated  into  the  solution.  The  solution  for  heat  sources  distributed  normal  to 
the  specimen  surfaces  is  obtained  by  using  the  discrete  plane  source  as  the  Green  function. 

Recently,  several  papers  have  treated  the  propagation  of  thermal  waves  in  linearly 
inhomogeneous  materials*’^  with  the  thermal  conductivity  (diffusivity)  gradient  normal  to  the 
specimen  surfaces.  From  a computational  point  of  view,  the  procedure  presented  in  this  note 
can  easily  treat  any  profile  as  a series  of  uniform  layers.  We  present  a sample  calculation 
for  a delaminated  diamond  film  on  a WC  substrate. 

The  procedure  is  also  applicable  to  solving  a wider  range  of  problems:  1)  Layered 
systems  containing  heating  sources  distributed  in  the  x-y  plane  parallel  to  the  layers.  The 
procedure  yields  the  x-y  spatial  Fourier  transform  of  the  temperature  distribution.  2) 

Layered  systems  having  a rectangular  cross  section.  The  procedure  yields  the  coefficients  of 
the  x-y  spatial  Fourier  series  expansion  of  the  temperature. 

2.  Matrix  formulation 

The  one  dimensional  homogeneous  differential  equation  that  describes  heat  diffusion 
for  a modulated  heating  source  of  angular  frequency  co  is. 


d^T  0) 

+i—T=0, 

dz^  ^ 


(1) 


where  D is  the  thermal  diffusivity  and  T is  the  complex  temperature  at  co.  The  time 
dependence  is  We  assume  that  D for  any  medium  in  the  problem  is  homogeneous, 

isotropic,  and  independent  of  T.  The  heat  source  has  not  been  included  in  the  equation; 
instead,  heating  is  treated  as  a boundary  condition.  The  solutions  can  be  written  as  the  sum 
of  exponentials  that  are  growing  (+)  and  decaying  (-)  in  the  positive  z direction. 
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r(z)=rv'«+r-^-'«, 


(2) 


where,  m^  = -i—,  and  and  T constants  to  be  determined.  The  solutions  to  eq.  (2)  can 
D 

also  be  expressed  in  the  form  of  waves,  called  thermal  waves,  if  we  make  the  substitution 
u=iw. 


The  temperature  can  be  treated  as  a two  dimensional  vector  T given  by 


T= 


T\z) 

[T-(z)) 


T'e- 


. Consider  a boundary  plane  at  between  media  a and  b,  as  seen 


in  figure  1,  with  no  heating  in  the  vicinity  of  the  boundary.  The  temperature  and  heat  flow 
are  continuous  across  the  boundary,  thus 


(3) 


(4) 


where,  is  the  limit  as  z approaches  the  boundary  from  the  right,  is  the  limit  as  z 
approaches  the  boundary  from  the  left,  and  kj  is  the  thermal  conductivity  of  medium  j.  The 
thermal  conductivity  is  related  to  the  thermal  diffusivity  by  K.=DpjCj,  where  pj  is  the  mass 
density  and  Cy  is  the  specific  heat.  Inserting  the  solution  eq.  (2)  into  eqs.  (3)  and  (4),  we 
obtain 


= r;(^)+r,(^),  (5) 

yX(^)  - y7;(^) = Y,r;(^)  - Y,r;(^) , (6) 


where  7y =«/«/• 

Equations  (5)  and  (6)  can  be  represented  as  a matrix  equation 

Ti=r,,T„  (7) 

where 


1 

2Y, 


Yi,-^Ya  YrY 

IY,-Y,  Y,+Y,J 


(8) 


Equation  (7)  relates  the  temperature  on  two  sides  of  an  interface. 

If  the  interface  between  two  media  contains  a thermal  resistance,  R [ref.  6],  then  eq. 
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(8)  becomes 


2Yi 


Yj-Y„-Yi,Y/ 
[.Yj-Y^+Y^Y^R  Y(,+Y„-Y(,Y/, 


The  relationship  between  the  temperatures  at  two  points 
separated  by  a distance  L=Zj-Zi  with  no  intervening  heat  source 


in  a given  medium  a 
is 


where. 


u«(^)= 


0 ^ 

0 


(8a) 


(9) 


(10) 


Thus,  eqs.  (7)  and  (9)  allow  us  to  express  the  temperamre  at  any  coordinate  as  a function  of 
the  temperamre  at  any  other  coordinate,  provided  there  is  no  intervening  heat  source. 

If  heat  is  applied  uniformly  to  the  boundary  at  a rate  ^ W cm"^,  then  eq.  (4)  becomes 


and  the  matrix  equation  becomes 


(11) 


Ti, 


=rz,A 


_g_ 

2Y, 


(12) 


If  the  heating  occurs  within  a single  medium,  eq.  (12)  simplifies  considerably  (r=l,  the 
identity  matrix);  thus,  when  performing  calculations,  it  is  most  convenient  to  consider 
heating  to  occur  within  a single  medium.  Heating  at  a boundary  between  two  different  media 
can  be  considered  as  heating  within  one  of  the  media  in  the  limit  that  the  distance  between 
the  heat  source  and  the  boundary  is  vanishingly  small. 

3.  One  dimensional  multilayer  heat  flow 

Consider  a multilayer  system  as  seen  in  figure  2.  We  want  to  solve  for  the 
temperamre  at  every  point  (or  plane)  for  heat  input  at  any  particular  point  (or  plane). 

Before  proceeding  with  the  solution,  let  us  define  some  nomenclamre.  If  we  have  N 
slabs  of  finite  thickness  with  media  to  the  right  and  to  the  left,  there  will  be  N+2  media  and 
N+1  interfaces.  We  begin  counting  media  from  the  left  starting  with  0 and  ending  with  N+\ 
to  the  right.  Each  interface  takes  the  label  of  the  medium  immediately  to  its  left.  The  origin 
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(z=0)  is  placed  at  interface  0.  The  coordinate  of  interface  j is  Zj  and  the  thickness  of  layer  j 
is  Lj.  The  nomenclature  for  the  temperature  vector  at  any  plane  z located  in  medium  j is 
where  f is  the  distance  from  the  left  boundary  of  medium  j,  that  is,  For 

brevity,  we  define  Tq=Tq{Q)  and  • 

Now,  let  us  consider  the  temperature  vector  relationship  in  the  vicinity  of  the  heat 
source.  Because  a single  medium  surrounds  the  source,  eq.  (12)  simplifies  to 


T/r)-T/r)= 


q 

2Y; 


(13) 


In  medium  0,  the  temperamre  must  remain  finite  in  the  limit,  z-»-oo,  thus,  Tq  (z)=0, 
whereas  in  medium  A^+1,  the  temperature  must  remain  finite  in  the  limit,  z~»oo;  thus, 

r^^i  = 0.  Hence  one  can  easily  show  that  and  are  both  proportional  to  Tq* 

and  Tj*{^+)  and  Tj~{^+)  are  both  proportional  to  If  Tj  \t-)=A  Tq  Tj  Tq  * 

and  then  eq.  (13)  becomes 


'■N+l 


B-ro"A=  — ^ 
2y. 


(14) 


/ N 

A^ 

/ \ 

B* 

where,  A= 

and  B = 

Equation  (14)  represents  two  simultaneous  equations  that 


can  be  solved  for  Tq  and 


The  solutions  are 


^ 2yjA^B--A-B^ 


(15) 


q A*+A 
2yj  A B -A  B^ 


(16) 


The  temperamres  in  all  of  the  layers  can  be  obtained  from  Tq*  and  ^iV+l  by  successive 
applications  of  the  vector  relationships  eqs.  (7)  and  (9). 


4.  General  Procedure 


To  solve  for  the  temperamre  distribution  throughout  the  multilayer  stack  for  a heat 
source  at  position  z=Zj_j+t,  we  must  first  solve  for  the  components  of  the  vectors  A and  B. 
These  are  given  by 
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(17) 


XTjj_j  X . . . Xr 5 2,1^^  1,0^ 

B =Uy(^-Ly)  X Tjj  ^j  X IJj  +j{-Lj  ^j)x  Tj  _^jj  ^2  X • • • 


• • • ^^N-2,N-l  ^N-l)  ^"^N.N+l  ^ 


(18) 


By  inserting  the  components  of  A and  B into  eqs.  (15)  and  (16),  we  obtain  the  temperatures 
at  the  outer  surfaces  of  the  multilayer  stack.  Successive  applications  of  eqs.  (7)  and  (9) 
allow  us  to  obtain  the  temperature  at  any  location  within  the  stack.  If  heating  occurs  at  an 
interface  between  two  media,  such  as  to  position  z=Zy_j,  we  just  take  the  limit  as 

If  there  are  heat  sources  at  two  or  more  coordinates,  we  solve  for  the  temperature  for 
each  source  separately  and  then  sum  the  solutions.  In  the  case  of  a heat  source  distributed 
along  z,  the  solution  obtained  above  can  be  used  as  the  Green  function  for  obtaining  the 
solution. 


5.  Spatially  Varying  Heat  Source  in  a Multilayer  System 

In  the  previous  discussion,  we  assumed  that  heat  was  applied  to  a plane  uniformly. 
Let  us  now  assume  a planar  heat  source  that  is  varying  in  the  plane  exists  at  z=f,  that  is, 
q=q{x,y).  The  heat  diffusion  equation  in  cartesian  coordinates  for  three  dimensional  heat 
flow  is  given  by 


q2t 

— + — + — +i— r=o. 

dx^  dy^  dz^  D 

The  general  solution  for  each  layer  in  a multilayer  system  is  of  the  form 


(19) 


Tj(.x,y,z)= f /{f;(k,,k^,z)+ 


(20) 


where. 

(21) 

f;k„k^,z)=cr(k,,ky‘‘\ 

(22) 

and 
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1/2 


(23) 


We  see  that  the  temperature  is  the  inverse  two  dimensional  Fourier  transform  of  two  terms. 
The  two  dimensional  Fourier  transform  of  q{x,y)  can  be  written  as 


(24) 


If  we  match  the  boundary  conditions  at  an  interface  between  two  layers  using  eqs.  (20)  and 
(24),  we  obtain  a matrix  equation  of  the  same  form  as  eq.  (12)  with  Fourier  transforms  of 
the  temperarnre  and  the  heat  source  replacing  the  temperature  and  heat  source  of  the  one 
dimensional  case.  In  addition,  eqs.  (21)  and  (22)  have  the  same  form  as  the  terms  in  eq.  (2). 
Thus,  the  formalism  developed  for  the  one  dimensional  case  can  be  seen  to  apply  directly  to 
the  Fourier  transforms  of  the  higher  dimensional  case.  In  the  case  of  an  axially  symmetric 
heat  source  in  the  xy  plane,  the  formalism  holds  as  well  except  that  we  must  use  Hankel 
transforms  instead  of  Fourier  transforms. 

6.  Multilayer  System  with  Rectangular  Cross-Section 

The  procedure  described  above  also  applies  to  a buried  planar  heat  source  in  a 
multilayer  system  having  a rectangular  cross-section  of  dimensions  and  L^.  In  this  case 
the  above  procedure  can  be  used  to  solve  for  the  coefficients  of  the  Fourier  series  analogous 
to  that  given  by  eq.  (15)  in  reference  4.  In  this  case  we  assume  heat  loss  from  the  edges  of 
the  films  is  neglible  compared  to  loss  through  the  end  faces  of  the  stack  and,  therefore,  can 
be  ignored. 

7.  Sample  calculation 

Thermal  waves  can  be  used  to  probe  for  delamination  of  a thin  film  on  a substrate 
and  the  one  dimensional  formalism  is  useful  for  analyzing  this  system.  The  system  consists 
of  an  air  ambient  (medium  0),  a film  (medium  1),  an  air  gap  (medium  2)  of  varying 
thickness,  and  a substrate  (medium  3);  N=2.  We  assume  that  the  modulated  heat  source  is 
in  medium  1 at  the  boundary  with  medium  0;  we  want  to  know  how  the  magnitude  and  phase 
of  the  temperature  at  the  film  surface  vary  as  a function  of  the  air  gap  thickness.  From  eqs. 
17  and  18  we  obtain 


2Yi  ’ 


(25) 


(26) 
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(27) 


(Yi-^Y2)(Y2-Y3)g'“^^(YrY2)(Y2n3)g“^  -.,t, 

4y,Y, 


and 


^-_(Yi-Y2)(Y2-Y3)^'“^'^(Yi~^Y2)(Y2~^Y3)g“^  u,l, 
4YiY2 


(28) 


Inserting  these  into  eq.  15,  we  obtain  the  surface  temperature  . 

The  magnitude  of  the  thermal  signal  is  given  by  | Tq""  [ , and  the  phase  is  given  by 

arctan(-Sro'/0iro').  Figure  3 shows  a plot  of  the  magnitude  and  phase  of  the  thermal  signal 
as  a function  of  air  gap  thickness  at  a modulation  frequency  of  200  Hz  for  a delaminated  film 
of  diamond  20^m  thick  on  a tungsten  carbide  substrate.  Table  1 lists  the  material  parameters 
used  in  the  computation. 

8.  Conclusion 

We  have  demonstrated  a straightforward  algorithm  for  solving  the  thermal  diffusion 
equation  in  one  dimension  for  a stratified  medium  with  a modulated  heating  source.  It  can 
be  applied  to  media  with  thermal  conductivity  gradients  if  the  gradient  is  normal  to  the 
specimen  surfaces.  In  this  case  the  media  can  be  decomposed  into  sublayers  of  uniform 
composition.  Interface  thermal  resistance  is  easily  incorporated  into  the  calculation.  If  the 
heat  source  has  a spatial  variation  in  the  x-y  plane  parallel  to  the  layers,  the  procedure  can  be 
used  to  obtain  the  x-y  spatial  Fourier  transforms  of  the  temperature  distribution.  In  the  case 
of  a mulitilayer  stack  of  rectangular  cross  section,  the  procedure  can  be  used  to  solve  for  the 
coefficients  of  the  x-y  Fourier  series  expansion  of  the  temperature.  A sample  calculation  for 
a delaminated  diamond  film  on  a WC  substrate  was  presented. 
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Table  1.  Numerical  Parameters  Used  in  the  Example 


medium 

/c  (W  cm*^  K'^) 

p (g  cm'^) 

C(Jg-'K-‘) 

0,  2 

0.00026 

0.00129 

1.01 

1 

9.6 

3.5 

0.51 

3 

0.95 

15. 

0.28 

z=^ 


h' 


Medium  a 


Medium  b 


Figure  1.  Heating  at  a boundary  between  two  media. 
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q 


Figure  2.  General  multilayer  system.  The  numbers  at  the  top  denote  the  boundary  indices, 
the  numbers  near  the  bottom  denote  the  media  indices,  Lj  denotes  the  layer  width,  and  q 
denotes  the  heat  input  a distance  f from  the  left  layer  boundary. 
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0 200  400  600  800  1000 
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Figure  3.  Magnitude  and  phase  of  the  thermal  signal  as  a function  of  air  gap  thickness  at  a 
modulation  frequency  of  200  Hz  for  a film  of  diamond  lOfim  thick  on  a tungsten  carbide 
substrate.  Table  1 lists  the  material  parameters  used  in  the  computation. 
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